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INTRODUCTION 

This work describes the application of unsteady 3-D Euler and 
Navier-Stokes equations to transonic flow past rotor blades, and 
wing-alone configurations. The computer code used in this study was 
developed under the U. S. Army Research Office support at Georgia 
Tech. The transonic wing-alone calculations were supported by the 
Lockheed Georgia Company under the I RAD program. 
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OBJECTIVES 


Methods based on the transonic small disturbance theory, and the 
full potential equation have matured to a point where they may be used 
by the industries for routine aeroelastic calculations. There is now 
the need to look at higher order techniques based on the Euler and the 
Navier-Stokes equations. The higher order solvers can serve two 
purposes. First, they provide a second estimate in situations where 
potential flow theory may fail (high transonic Mach numbers, strong 
shock waves), and provide benchmark runs for the validation of 
potential flow codes. Secondly, they allow the designer to study 
phenomena such as high angle of attack transonic maneuvers, supersonic 
fighter aerodynamics, and 3-D separated flow around highly loaded 
rotor blades. 


1. To Describe a Solution Procedure for 
the Numerical Solution of the 3-D 
Compressible Viscous or Inviscid Flow 

2. Apply this procedure to a number of 
fixed and rotary wing problems of 
interest 


Figure 1 



GOVERNING EQUATIONS 


The equations governing three-dimensional unsteady compressible 
flow are the Navier-Stokes equations. If the viscous terms are 
neglected, the Euler equations result. The present solution techniques 
are designed to work efficiently with both the Navier-Stokes and the 
Euler equations. All calculations have been done on an algebraically 
generated body-fitted coordinate system (c,n»C) * which is allowed to 
move with time and follow the motion of the solid. The flow 
properties of interest at a given time level are p : the density, 
u,v,w: the velocity of the fluid in an inertial coordinate system, and 
e: the total energy of the fluid per unit volume. The quantities U,"V 
and W are the contravariant components of velocity along the 5- , n~ 
and £- directions respectively. Also, p is the pressure. 
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GOVERNING EQUATIONS (CONTD.) 


The quantities c x » C y > 5 Z etc * are the metrics of the 
transformation computed numerically using standard second order 
accurate finite difference formulas. The quantity J is the Jacobian 
of transformation. 


NAVIER-STOKES EQUATIONS 
q r + (E-E v ) { + (F-F v ), + (G-G v )j = 0 
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TURBULENCE MODEL 


A two-layer eddy viscosity model developed by Baldwin and Lomax 
is used in this study (Ref. 1). For the mildly separated flows 
considered here this model has proved adequate. Here 1 is the mixing 
length in the inner layer, proportional to the distance m from the wall, 
and the van Driest damping factor. In the outer layer, F max is a 
measure of the velocity scales within the shear layer, while d max 
is a measure of the length scale. At large distances from the shear 
layer the eddy viscosity is designed to approach zero through the use 
of the i ntermi ttency factor F^. 


• Baldwin-Lomax two-layer algebraic model used for 
eddy viscosity. 

• Inner Layer: 

m T = for d < d c 

= (Kd) t 1 “ exp(-d + /A + ) ] 

d + = d(pr w ) ,/2 A< 

• Outer Layer: 

= '.OlGSpc^F^ for d > d c 

F w “ min(d max F max- C 2 d mx0 2/F max > 
F max = max tV' /K ] 

F k = [1 + 5 - 5 (C3 d/d max )6rl 
Q = max|v| - min|v| 


Figure 4 
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HYBRID ALGORITHM 


In the present procedure, the time derivatives appearing in the 
governing equation are discretized using a two-point backward 
difference formula. The derivatives along the £- and £ directions have 
been kept at the new time level (n+1) where the solution is sought. 
The spanwise (n~) derivatives have been explicitly evaluated using the 
latest available information at the inboard station during odd time 
steps and outboard station during the even time step. Thus, the 
computational stencil resembles a plane Gauss-Seidel algorithm, where 
the spanwise sweeps are performed in opposite directions on alternate 
iteration levels. It may be shown from a linear stability analysis 
that this technique leads to a stable algorithm. 


Implicit Euler rule: 


r 1 - q" ♦ At|? 


n+1 


Evaluate spanwise term explicitly: 

q n+1 = q 11 - At(8|E n+1 + 6 JJ F n,n+1 + 8^G n+1 ) 

where spanwise term alternates between: 
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TIME LINEARIZATION 


The fact that some of the quantities are to be computed at the 
new time level (n+1) means that the resulting system of equations is 
algebraic, but highly nonlinear. To avoid an iterative solution of 
non-linear equations, the well known Beam-Warming linearization (Ref. 
2) is applied to the flux terms along the and c directions. The 
result is a system of linear, block pentadiagonal equations in which 
the unknown is the 'delta 1 change in the flow properties between 
adjacent time levels. 


Second-order expansion: 

E n+1 = E n + [A n ](q n+1 - q n ) + 0(At 2 ) 

G n+1 = G n + [C n ](q n+1 - q n ) + 0(At 2 ) 

where A and C are the Jacobian matrices: 

[A] = 3E/3q and [C] = 3G/3q 

The following linear system results: 

[I + At ( 6 ^A n + 6^C n ) ]Aq Il+1 = R n ' n+1 

Aq n+1 * q 11 * 1 - q* 

R n,n+1 = - At( 5^ (E n - E J» n+1 )+6 ^ (F -F v ) n,n+1 +6 s (G n -G^ n+1 )) 


Figure 6 
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APPROXIMATE FACTORIZATION 


The direct inversion of the pentadiagonal block matrix equations 
is expensive. In literature, a number of techniques are available, 
based on the approximate factorization techniques such as L-U 
decomposition or ADI decomposition. The purpose of these techniques is 
to break up the coefficient matrix into smaller, easily inverted 
matrices. In this work an ADI factorization is used to arrive at a 
system of block tri diagonal equations, which may be inverted using the 
Thomas algorithm. 


ADI Solution in the Airfoil Plane 


Approximate factorization: 

[I + At(6 f A + 6^C) ] * [I + AtS^A] [I + At6jC] + 0(At 2 ) 

Gives two linear systems with block tridiagonal matrices: 

[I ♦ At6 { A]Aq* n+1 = R n ' n+1 
[I + At6^C]Aq n+1 = Aq* n+1 

At every time step, the solver marches through the radial 
stations explicitly performing two matrix-inversion sweeps at 
each station. 


Figure 7 
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ARTIFICIAL DISSIPATION TERMS 


The use of pure central differences to advance hyperbolic, or 
weakly parabolic equations can lead to numerical instabilities. A 
variety of artificial dissipation forms have been suggested in 
literature to overcome this instability. In this work, a fourth order 
explicit dissipation form is used, and to allow large amounts of 
explicit dissipation to be used without instability, a second order 
implicit dissipation term is added to the left side. 


• Using central differencing alone leads to odd-even 
decoupling. 

• Nonlinearities of the equations produces high- 
frequency errors which grow. 

• Second-order implicit dissipation and fourth-order 
explicit dissipation used. 

• Fourth-order implicit dissipation stabilizes the 
scheme even more, but results in penta-diagonal 
systems. 

• It can be shown that the dissipation results in an 
upwinded scheme. 


Figure 8 
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FINAL FORM OF THE DIFFERENCE EQUATIONS 

The final form of the discretized equations, used in the computer 
code are shown below. 


[I + At6gA - Ate J -1 V^A^J]Aq* n+1 = R n ' n+1 _ D fi n ' n+1 
[I + AtS^C - At ( e ( +e^) J** 1 V^A^j]Aq n+1 = Aq* n+1 

where, the explicit dissipation is: 

D e n ' n+ 1 = Ate E J _1 [(V^) 2 + (V^) 2 + (V s A s ) 2 ]jq n 

and the variable implicit coefficient is: 

+ «!><**+ 5y + £ 

As an example, 

<V £ A { ) 2 Jq" ■= (q it2 - 4q i+1 ♦ 6q A - 4^ + 

(V { A { )Jq n = (q i+1 - 2q, ♦ « i . 1 ) n 

Figure 9 
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TREATMENT OF ROTOR WAKE 


In helicopter applications, the treatment of the shed tip 
vortices requires special considerations. The finite difference grid 
is usually not large enough or fine enough to capture the many 
revolutions of the tip vortices shed by several blades. In this work, 
only the portion of the shed vorticity immediately downstream of the 
rotor blade is captured by the finite difference scheme. The rest of 
the vorticity, and the wake behind the other blades is kept track of 
using a Lagrangean approach. This approach was first proposed by Tung 
and Caradonna (Ref. 3). 


• Downwash due to tip vortex significantly affects the lift 
of the blade 

• Resolving the tip vortex by finite difference techniques 
is not possible with current computer resources 

• Effects of tip vortices lying outside of the computational 
domain must be included 



• Use CAMRAD or some other wake code to obtain 
effective partial angle of attack distribution 


Figure 10 



LIFTING VISCOUS FLOW PAST A HOVERING ROTOR 


As a first application of the solution technique described, the 
subsonic lifting flow past a two bladed rotor system in hover is 
considered. The blades were made of NACA 0012 airfoils and had a 
rectangular planform. The collective pitch was 8 degrees, and the tip 
Mach number was 0.44. There is an extensive set of experimental data 
available for this configuration (Ref. 4). Here, the computed pressure 
distrbution at a number of radial stations is plotted and compared 
with experimental data. Good agreement is observed at all locations. 



Figure 11 
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Lift Coefident 


SPANWISE LOADING ON A HOVERING ROTOR 

The integrated lift distributions are plotted along the rotor 
radius below. It is seen that good agreement with experiments is 
found. 

These calculations were done on a 79 x 23 x 45 grid with 50 
points at each radial location on the rotor, 11 radial locations on 
the rotor and 45 points in the normal direction. They required 3.9 
seconds per time step on the CRAY XMP. 


Lift Coefficient Distribution 



Figure 12 
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NON-LIFTING TRANSONIC FLOW OVER AN ADVANCING ROTOR 


As a second example, the transonic viscous flow past a rotor 
blade tested at ONERA is considered. The tip Mach number in this case 
was 0.6, while the advance ratio (forward speed/tip speed) was 0.45. 
In the next several pages, comparisons between experiments and the 
Navier-Stokes solutions are given at the 84% span station. For the 
sake of completeness some Euler solutions are also shown. It is seen 
that both the Euler and the Navier-Stokes solutions give acceptable 
agreement with experiments. It is also seen that the Navier-Stokes 
results predict the shock locations and strength more accurately. 
These calculations were done on a 121 x 19 x 45 grid and required 4.6 
sec per iteration on a CRAY XMP. Several thousand time steps were 
needed to advance the solution from zero degree azimuth to 360 degree 
azimuth in time, (through a full revolution). 




Figure 13 
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LIFTING TRANSONIC FLOW PAST A RECTANGULAR WING 


As a final application of the present approach, the unsteady 
transonic flow past a rectangular supercritical wing tested at NASA 
Langley Research Center is presented. The freestream Mach number was 
0.7 and the mean angle of attack was 1.98 degrees. The wing was 
constrained to oscillate in pitch at a frequency of 10 Hz. In the 
following figures, the in-phase and the out-of-phase components of the 
surface pressure distribution are plotted at several locations. 
Overall, a reasonably good agreement is observed. 



FIGURE 14. R376E21 INVISCID REAL AMPL=1.044 2000 S/C 5000 S , M = 
0.70, a 0 = 2.0 deg. 


371 


- 10.0 0.0 10.0 - 10.0 0.0 10.0 


LIFTING TRANSONIC FLOW PAST A RECTANGULAR WING (concluded) 
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FIGURE 14 . R376E21 INVISCID IMAGINARY AMPL=1.044 2000 S/C 5000 S 

(Concluded) 
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CONCLUDING REMARKS 


A promising approach for the numerical solution of 
three-dimensional Euler and Navier-Stokes equations has been 
described. Additional work is needed for improving the efficiency of 
the present procedure. It is hoped that the techniques presented here 
will find use in fixed and rotary wing aircraft analysis. For 
additional studies and code correlations the reader is referred to 
Refs. 5-8. 


1. A solution technique for the 3-D 
Euler and Navier-Stokes equations 
has been developed. 

2. A number of interesting fixed and 
rotary wing applications have been 
presented. 

3. Additional work towards improving the 
solution efficiency is now underway. 


Figure 15 
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